/***************************************************************
                        cb_compute_order.c 
Routines that compute the order of which clones are added to
the cbmap

Zsel_zap - user selects clones from contig to assemble
Zcalc_zap - assemble clones in contig
Zagain    - reassemble clones in CBmap
Zbuild_contigs__zap - assembled contigs with all clones

Do_BestOf - used by these routines to try N solutions based on diff clones
Greedy_make_Zset - finds the next set (CBmap) of clones to use. 

Changes to results since v7.2: 
1. When Ntries goes through all unburied and then starts with buried,
it now goes through all buried also (or until Ntries); it used to just keep 
repeating the last one over and over. The makes average and low different
2. The Avg Score at the end was the average of the average scores instead of
the best scores.
3. A change in cb_assemble.c to do with buried clones will make the Qs and Score
slightly different.

DO_SHARED_NTRIES = (num_clones*Pz.bestInt) > PARALLEL_CUTOFF)
***************************************************************/
#include <stdio.h>
#include <malloc.h>
#include <unistd.h>
#include <stdlib.h>
#include <assert.h>
#include <semaphore.h>

#include "clam.h"
#include <sys/types.h>
#include <sys/wait.h>
#include <sys/times.h>

extern int Num_threads;
extern int total_cmp, total_cnt;
extern int cnt_overlap, cnt_nodes;
extern void clearMsg(), unbury();
extern void quitCB(), showClam(), showClam2(), settime(), scanClamB();
extern void cleanCB(), cbmap_sort(), OK_Zset(),  Zproj_results();
extern int Zbuild_matrix();
extern void Zmark_possible_buried(), ZbuildCpM(), Zfree_cb_data_structures();
extern ZPQZ *insertZpq();
extern int init_Zclone(), get_Zset();
extern void free_Zclone();
extern void Zprepare_for_display(), init_Zset(), make_this_Zset();
extern void big_wrap_up(), little_wrap_up();
extern int shared_node_creation();
extern void pr_end_timer();

int  Do_BestOf();
int Greedy_make_Zset();
int cnt_cmp;
int parallelize_ntries(int best, float *save_score, int *save_j,
                     float * avg_high_score, int * save_c);
int get_seed_clone(int);
static int assign_clones_to_CBmaps();
static int cnt_par;

static void init_simulation();
static int compute_chimerics();
static void print_simulation_results();
static int compute_out_of_order();

/***************************************************************
         Def: Zbuild_contigs_zap()
****************************************************************/
void Zbuild_contigs_zap()
{
int i, c, save, ctg;
CLONE *clp;
clock_t start=0,end=0, start2=0, end2=0;
struct tms tmsstart,tmssend, tmsstart2, tmssend2;

  init_simulation();
  
  if (graphActivate(g99)) quitCB();
  else cleanCB();
  if (contigs[0].count==0) {
     sprintf(ZBuf,"No singletons"); showClam2();
     return;
  }
  clearMsg();
  LastCB=Build;
  strcpy(Fn,"Build");

  for (i= contigs[0].start; i != -1; i = clp->next)
  {
     clp = arrp(acedata, i, CLONE);
     if (clp->clone[0] == '!') continue;
     if (clp->match[0]!=' ') unbury(i);
  }

/**********BUILD SPARSE MATRIX ***************/
  if (!Zbuild_matrix(0)) {
      printf("FPC error: cannot build sparse matrix\n");
      return;  
  }
  
  start=times(&tmsstart);
  if(SHARED_OPTION && DO_SHARED_NTRIES(contigs[0].count))
  {
     if(shared_node_creation() < 0) {
          fprintf(stderr, "FPC ERROR: shared node creation\n");
          goto bail;
     }
     printf("\nComplete building parallel sparse matrix...\n\n");
     fflush(0);
  }
  else {
     for (i=0; i<ZZ.size; i++)
     {
         if (!Zcreate_node(i)) {
            fprintf(stderr,"FPC ERROR: serial node creation\n");
            goto bail;
         }
     }
     printf("\nComplete building serial sparse matrix...\n\n");
     fflush(0);
  }
  end = times(&tmssend);
        
  Zmark_possible_buried();

/********* GREEDY ORDER CLONES ****************/
  cnt_par=0;
  start2=times(&tmsstart2);
  if (assign_clones_to_CBmaps()==0) return;

  save = currentctg = ctg = getcontigsmax()+1;

  for (ZS=0; ZS<n_Zset; ZS++) {
                                   /* Do_BestOf calls parallel */
     Do_BestOf();
     OK_Zset(0, ctg);
           /* this is only necessary to free memory **/
     for (c=0; c< Zset[ZS].n_clone; c++) {
         i = Zset[ZS].csort[c].i;
         free_Zclone(i);
         arrp(acedata, ZZ.matrix[i].cin, CLONE)->oldctg = ctg;
     }
     currentctg = ctg = getcontigsmax()+1;
  }

/************ PRINT STATS and CLEAN UP *****************/
bail: 
  end2 = times(&tmssend2);
  currentctg=0;
  ZBuf[0]='\n'; 
  big_wrap_up(0);

  if (n_Zset>0) {
    int ctgcnt[4], max=0;
    for (i=0; i<4; i++) ctgcnt[i]=0;
    for (i=0; i< n_Zset; i++) {
        max = MaX(max, Zset[i].n_clone);
        if (Zset[i].n_clone > 50) ctgcnt[0]++;
        else if (Zset[i].n_clone > 25) ctgcnt[1]++;
        else if (Zset[i].n_clone > 3) ctgcnt[2]++;
        else ctgcnt[3]++;
    }
    sprintf(ZBuf,
   "Create %d contigs (%d:%d): Max %d, %d (>50), %d (50:26), %d (25:4), %d (3:2)\n",
      n_Zset, save, save+n_Zset-1,
      max, ctgcnt[0], ctgcnt[1], ctgcnt[2], ctgcnt[3]); Outlog;
    if (cnt_par > 0) printf("Run layout in parallel for %d contgs (cutoff %d)\n",
           cnt_par, PARALLEL_CUTOFF);
    pr_end_timer("NxN Pairs:", end-start,&tmsstart,&tmssend);
    pr_end_timer("Layout:   ", end2-start2,&tmsstart2,&tmssend2);
    printf("\n");
    sprintf(ZBuf,"Create %d contigs",n_Zset);
  }
  else strcpy(ZBuf,"No contigs.");
  showClam2();

  if (SIM_FLAG)
       print_simulation_results();

  settime(&Bd.date);
  ZbuildCpM(); /* set table values for saving in fpc file */
  if (!Zbatch_flag) ClamCtgDisplay2();
  Zfree_cb_data_structures();
}
/****************************************************************
                       Def: Zcalc_zap
This is called by Calc in the Contig Analysis (clam.c)
And by the DQer and Reanalysis in the Main Analysis (c2am.c)
And Ends->Ends (ends2ends.c)
****************************************************************/
void Zcalc_zap()
{
int i, save;
   save = LastCB;
   scanClamB(); // otherwise you have to hit return WMN
   cleanCB();
   if (save==DQ) strcpy(Fn,"DQer");
   else if (save==NoQ && Pz.rebuild) strcpy(Fn,"Q-");
   else if (save==NoQ && !Pz.rebuild) strcpy(Fn,"Q~");
   else if (save==EndMerge) strcpy(Fn,"End");
   else if (save==IBC) strcpy(Fn,"IBC");
   else strcpy(Fn,"Calc");
   LastCB=Calc;  /* LastCB was set in DQ and Q- so to know how to write Fn */

   if (!qstuff) printf("\n");
   Zadd = Zorder = 0;

/************ BUILD SPARSE MATRIX *******************/
   if (!Zbuild_matrix(currentctg)) return;
   
   if(SHARED_OPTION && DO_SHARED_NTRIES(contigs[currentctg].count))
   {
      if(shared_node_creation() < 0) {
         fprintf(stderr, "FPC ERROR: shared node creation\n");
         return;
      }
   }
   else 
   {
      for (i=0; i<ZZ.size; i++) 
         if (!Zcreate_node(i)) return;
   }
   Zmark_possible_buried();

/*********** GREEDY COMPUTE ORDER ******************/
  if (assign_clones_to_CBmaps()==0) return;
  for (ZS=0; ZS<n_Zset; ZS++) {
       Do_BestOf();
  }
  big_wrap_up(1);
  LastCB = save;
  if (!qstuff) printf("Complete Calc (use OK to instantiate)\n");
}
/****************************************************************
                       Def: Zsel_zap
The user select a set from the contig to be assembled.
25 mar 05 - this used to set cutoff to 1e-01 so all would go into one contig.
But can cause hanging if done on large selected set (>1000).
So user will have to lower cutoff it they want this behavior.
****************************************************************/
void Zsel_zap()
{
int i, f;
CLONE *clp;

   for (f=0, i= contigs[currentctg].start; i != -1; i = clp->next) {
     clp = arrp(acedata, i, CLONE);
     if (clp->selected==1) f++;
   }
   if (f<=1) {
      sprintf(ZBuf,"** Select 2 or more clones."); 
      showClam(); 
      Outlog;
      return;
   }
   printf("Selected %d clones\n", f);
   
   cleanCB();
   strcpy(Fn,"Select");
   LastCB=Select;
   Zadd=Zorder=0;
   qstuff=0;

/************ BUILD SPARSE MATRIX *******************/
   if (!Zbuild_matrix(currentctg)) return;
   
   if(SHARED_OPTION && DO_SHARED_NTRIES(contigs[currentctg].count))
   {
      if(shared_node_creation() < 0) {
         fprintf(stderr, "FPC ERROR: shared node creation\n");
         return;
      }
   }
   else 
   {
      for (i=0; i<ZZ.size; i++) 
         if (!Zcreate_node(i)) return;
   }
   Zmark_possible_buried();

/************* Order clones ******************/
   if (assign_clones_to_CBmaps()==0) return;
   for (ZS=0; ZS<n_Zset; ZS++) Do_BestOf();
   big_wrap_up(1);
}
/****************************************************************
                       Def: Zagain
****************************************************************/
void Zagain_zap()
{
int i, save, cnt=0;

   sprintf(Fn,"Again");
   save = Zadd;
   Zadd=Zorder=0;
   
   for (i=0; i<ZZ.size; i++) {
       if (ZZ.matrix[i].cbmap == ZS) {
           if (ZZ.matrix[i].grand <= 3) cnt++;
           init_Zclone(i);
       }
   }
   if (cnt >= Zset[ZS].n_clone-1) {
      for (i=0; i<ZZ.size; i++) 
           ZZ.matrix[i].grand = ZZ.matrix[i].save_grand;
   }
   Greedy_make_Zset(get_seed_clone(ZS));  /* never parallel */
   if (n_Zset>0) cbmap_sort();
   Zadd=save;
   little_wrap_up();
}

/******************************************************************
                DEF: GREEDY_MAKE_ZSET()
*****************************************************************/
int Greedy_make_Zset(int C)
{
    ZPQZ *pq;

   if (C == -1) return 0;
   
   init_Zset(ZS, 1);
   insertZpq(C, 1000);  
   while ((pq = popZpq())) 
   {   
       make_this_Zset(pq);  /* add to CBmap */
   }
   if (Zset[ZS].n_clone <= 1) { /* shouldn't happen, but not fatal */
       fprintf(stderr,"FPC warning - one clone contig - ignored\n");
       return 0;
   }
   Zprepare_for_display();
   if (!Zbatch_flag && !SHARED_OPTION && graphInterruptCalled()) {
          fprintf(stderr,"F4 pre-mature termination of CALC or Build Contigs\n");
          n_Zset=0;
          return 0;
   } 
   return 1;
}

/*****************************************************************
                DEF: DO_BESTOF
called by IBC, Calc, and Build Contigs

Greedy_this_Zset is called once before calling here in order
to initialize ZS and Zet[ZS] 
****************************************************************/
int Do_BestOf() 
{
int j, c, s;
int  save_c=0, save_j=0, nTries=1;
float save_score=0.0 , avg_high_score =-1;
char p=' ';

   if (Zset[ZS].n_clone>1000) {
       printf("Building large cbmap of %d clones....\n",Zset[ZS].n_clone);
       fflush(0);
   }

   avg_high_score = Zset[ZS].avg = Zset[ZS].score = 0.0; 
   Zset[ZS].low = 99999;

   nTries = MiN(Pz.bestInt, Zset[ZS].n_clone-1);
  
   if(Num_threads <=1 || nTries <=2 || !DO_SHARED_NTRIES(Zset[ZS].n_clone)) 
   { 
      for (j=0; j<nTries ; j++) { /* next n tries */
           for (c=0; c<Zset[ZS].n_clone; c++) 
                init_Zclone(Zset[ZS].csort[c].i);

           if ((s = get_seed_clone(ZS)) == -1) break;
           if (!Greedy_make_Zset(s)) {
               fprintf(stderr,"FPC error in greedy algorithm\n");
               goto done;
           } 
           if (Zset[ZS].score > save_score) {
                  save_score = Zset[ZS].score;
                  save_c = Zset[ZS].csort[0].i;
                  save_j = j;
           }
           Zset[ZS].avg += Zset[ZS].score;
           if (Zset[ZS].low>Zset[ZS].score) Zset[ZS].low=Zset[ZS].score;
      }
   }
   else
   {
       if((parallelize_ntries(nTries, & save_score, & save_j,
          & avg_high_score, & save_c))==0) 
               return 0;
       Zset[ZS].avg = avg_high_score;   
       p = 'p';
       save_j = -1; /* not returning corect value */
       cnt_par++;
   }
   
   if (save_j != nTries-1) { 
        for (c=0; c<Zset[ZS].n_clone; c++) 
            init_Zclone(Zset[ZS].csort[c].i);
        if (!Greedy_make_Zset(save_c)) return 0;
        if (save_score != Zset[ZS].score) 
           printf("FPC warning: last score %6.3f != current %6.3f\n",
                  save_score, Zset[ZS].score);
    }
    Zset[ZS].avg /= (float) nTries;

done: ;
         /* the mult and div is so give the same round off as will be stored */
   if (!SIM_FLAG) {
    sprintf(ZBuf,
" %c%s Ctg%-4d CBmap %-4d Clones %-4d CBs %-5d  Qs %-3d Score %4.3f Avg %4.3f Low %4.3f\n",
          p, Fn, currentctg, ZS+1, Zset[ZS].n_clone, Zset[ZS].n_band, Zset[ZS].n_Qs,
          ((Zset[ZS].score*1000.0)/1000.0), ((Zset[ZS].avg*1000.0)/1000.0), 
          ((Zset[ZS].low*1000.0)/1000.0)); Outlog;
   }
    return 1;
}

/************************************************************
   Def: parallelize_ntries
  Called from Do_BestOF : Used for parallelizing "n tries"
*************************************************************/

int parallelize_ntries(int nTries, float *save_score, int *save_j, 
                      float * avg_high_score, int * save_c)
{
    int i,j,index,status, best_thread = -1; 
    int iter_per_thread, extra_iter_thread; 
    int exit_status = -1;
    int pid=0,thread_number,save_iteration,c;
    float high_save_score = -1;
    struct Ntries_thread * Nallthread;
    int *seedmem, offset=0;
#define NThreads 50
    int fd[NThreads][2];
    int *seeds[NThreads];
    
    if (NThreads < Num_threads) {
       fprintf(stderr,"Only allowed %d threads. Change in cb_compute_order\n",
              NThreads);
       return 0;
    }
    *avg_high_score = Zset[ZS].score;
   
    Nallthread=(struct Ntries_thread *) malloc (sizeof(struct Ntries_thread)
                                             * (Num_threads)); 
    NOMEM3(Nallthread, "Nallthread", 0);
    
    if(nTries <= Num_threads) nTries = Num_threads;
    iter_per_thread = nTries / Num_threads;
    extra_iter_thread = nTries % Num_threads; 
     
    /* Fill up thread data structre for number of tries/thread */
    for(i=0;i<Num_threads;i++)
    {
       Nallthread[i].ntries = iter_per_thread;
       Nallthread[i].save_score=0;
       Nallthread[i].avg_score=0;
       Nallthread[i].save_j=0;
       Nallthread[i].save_c=0;
       Nallthread[i].low= 9999.9; 
       Nallthread[i].pid= 0;  /* WMN prevent valgrind complaining about uninitialized var */
       if(extra_iter_thread >0)
       {
          Nallthread[i].ntries ++;
          extra_iter_thread --;
       }
    }   

       /* for intialization of the starting seeds*/
    seedmem =(int*) (malloc(sizeof(int)*nTries+1));
    NOMEM3(seedmem, "seedmem", 0);
    for(i=0;i<Num_threads;i++)
    {
        seeds[i] = &seedmem[offset];
        for(index=0; index < Nallthread[i].ntries;index++) 
           seeds[i][index] = get_seed_clone(ZS); 
        offset += Nallthread[i].ntries;
    }
 
    i=0;

    /* Creating the children and the associated pipes */
    while(i<Num_threads)
    {
        if (pipe(fd[i])<0)
        {
            fprintf(stderr, "FPC Error creating pipe\n");
            return 0;
        }
        if((pid=fork())>0) /* parent */
        {
           Nallthread[i].pid=pid;
           i++;
        }
        else if(pid<0) {
           fprintf(stderr,"\n\n ERROR....cannot create threads.....\n\n");
           return 0;
        }
        else               /* child */
          break;
    }

/***********  CHILD CODE ***********/
    if (pid==0)
    {
       thread_number=i; 
       close (fd[i][0]);

       for (j=0;j<Nallthread[thread_number].ntries;j++)
       {
          if (Zset[ZS].n_clone>1000) {
             printf("Thread %d try %d/%d           \r",thread_number,j+1,Nallthread[thread_number].ntries);fflush(0); //WMN
          }
          save_iteration = j * Num_threads + thread_number + 1;

          for (c=0; c<Zset[ZS].n_clone; c++)
             init_Zclone(Zset[ZS].csort[c].i);
 
          if (!Greedy_make_Zset(seeds[thread_number][j]))
          {
              fprintf(stderr,"FPC error on thread %d during Greedy\n",thread_number);
              write(fd[thread_number][1],&Nallthread[thread_number],
                   sizeof (struct Ntries_thread)); 
              close (fd[thread_number][1]);
              _exit(2);
          }

          if (Zset[ZS].score > Nallthread[thread_number].save_score) {
              Nallthread[thread_number].save_score = Zset[ZS].score;
              Nallthread[thread_number].save_c = Zset[ZS].csort[0].i;
              Nallthread[thread_number].save_j = save_iteration;
          }
          Nallthread[thread_number].avg_score += Zset[ZS].score;
            
          if (Nallthread[thread_number].low>Zset[ZS].score) 
          {
             Zset[ZS].low = Zset[ZS].score;
             Nallthread[thread_number].low = Zset[ZS].low; 
          }
      } 
      if (Zset[ZS].n_clone>1000) {
         printf("Thread %d done           \n",thread_number); //WMN
         fflush(0);
      }     

      /* Send data back to the parent */
      write(fd[thread_number][1],&Nallthread[thread_number],
               sizeof (struct Ntries_thread)); 
      close (fd[thread_number][1]);
       _exit(0);  
    }  /*end for child thread*/
   
/*********** PARENT CODE ***********/

    /* Close write ends of pipe for parent */
    for(i = 0; i< Num_threads; i ++)
    {
        close (fd[i][1]);
    }

    /*Make the parent wait till all children complete*/
    i=0;   
    while(i++<Num_threads)
    {
        exit_status=-1;
        pid=wait(&status);

        if(WIFEXITED(status))
           exit_status=WEXITSTATUS(status); 

        if (exit_status==2)
        {
           for (j=0;j<Num_threads;j++)
           {
              if (Nallthread[j].pid!=pid)
                 kill(Nallthread[j].pid,0);
           }
           printf("\n returning due to premature\n");
           free(Nallthread);
           free(seedmem);
           return(0);
        }

        /* find pid's index into Nallthread */
        for( j=0; j < Num_threads ;j ++)
        {
            if (Nallthread[j].pid == pid)
                break;
        }

       /*read data from the child */
       read(fd[j][0],(char *)&Nallthread[j],
                  sizeof (struct Ntries_thread));
       close(fd[j][0]); 

       /* the child did not have pid set, so lost this value on above read */ 
       Nallthread[j].pid=pid;
    }
    /** All threads are done, find avg, low and high **/
    for(i=0;i<Num_threads;i++)
    {
        *avg_high_score+= Nallthread[i].avg_score;     

        if(Nallthread[i].save_score > high_save_score )
        {
            high_save_score = Nallthread[i].save_score;
            best_thread = i;
        }
        if (Zset[ZS].low > Nallthread[i].low)
            Zset[ZS].low = Nallthread[i].low; 
    }    

    *save_j = Nallthread[best_thread].save_j;
    *save_c = Nallthread[best_thread].save_c; 
    *save_score = Nallthread[best_thread].save_score;
 
    free(Nallthread);
    free(seedmem);
    return(1);
}

/******************************************************
                 Def: get_seed_clone
 Called from parallelize_ntries. Used to fill up 
 the clones to start with data structure for every thread
********************************************************/
int get_seed_clone(int zset)
{ 
    int C,grand,c,nb,i;
          
    /* find clone to start with */
    /* buried have grand of 1, so are only used if all others used in tries */
    for (C=-1, grand=1, nb=c=0; c< ZZ.size; c++)    
    {
        if (ZZ.matrix[c].cbmap == zset)
        {
             i=0;  
             if (ZZ.matrix[c].grand >= grand) i=1;
             else if (ZZ.matrix[c].grand == grand) {
                  if (ZZ.matrix[c].nbands > nb) i=1;
                     /* check the name so the solution is 
                        the same regardless of order of input clones */
                  else if (ZZ.matrix[c].nbands == nb && 
                         strcmp(arr(acedata, ZZ.matrix[c].cin, CLONE).clone,
                         arr(acedata, ZZ.matrix[C].cin, CLONE).clone)>0)  i=1;
             }
             if (i) {
                   nb = ZZ.matrix[c].nbands;
                   grand = ZZ.matrix[c].grand;
                   C = c;
             }
        }
    }
    if (C == -1)  return -1; 
//printf("%d %d %d\n", zset, C, ZZ.matrix[C].grand);
    ZZ.matrix[C].grand=0;
    return C;
}
/***************************************************************************
                       DEF: assign clones to CBmaps
***************************************************************************/
static int assign_clones_to_CBmaps()
{
int z, c, i;
ZPQZ *pq;
void add_nodes_to_pq();

         /* all clones in ZZ.matrix are assigned a cbmap */
   while ((c = get_seed_clone(-1)) != -1) {
        if ((ZS=get_Zset()) == -1) {
           fprintf(stderr,"FPC error: cannot get zset\n");
           return 0;
        }
        ZZ.matrix[c].cbmap = ZS;
        ZZ.matrix[c].pick = 1;   /* will not be reinserted into Q */
        add_nodes_to_pq(c);
        Zset[ZS].n_clone++;
        while ((pq = popZpq())) 
        {   
            ZZ.matrix[pq->C].cbmap = ZS;
            ZZ.matrix[pq->C].pick = 1;
            add_nodes_to_pq(pq->C);
            Zset[ZS].n_clone++;
        }
   }
        /* the Zset structure stores the clone list in csort */
   for (z=0; z<n_Zset; z++) {
        Zset[z].csort = (struct Zset_csort *) 
                 malloc(Zset[z].n_clone * sizeof(struct Zset_csort));
        NOMEM3(Zset[z].csort, "struct Zset_csort",0);
        for (i=c=0; i<ZZ.size && c<Zset[z].n_clone; i++)  
            if (ZZ.matrix[i].cbmap == z) { 
                ZZ.matrix[i].grand = ZZ.matrix[i].save_grand; /* reset for greedy */
                Zset[z].csort[c].i = i;
                c++;
                if (Zset[z].n_clone < c) 
                   printf("add_clones_to_cbmap: %d %d\n",z,c);
            }
   }
return 1;
}




//This function prints the results of the simulation

static int count_singletons()
{
    int i;
    int sings = 0;
    CLONE* clp;
    for(i=0;i<arrayMax(acedata); i++)
    { 
        clp =  arrp(acedata, i, CLONE);
        if (clp->clone[0] == '!') continue;
        if(clp->ctg == 0 ) sings++;
    }
    return sings;
}
static void print_simulation_results()
{
    int i =0;
    int num_Q = 0, num_Q_ctg=0;
    float avg_score = 0;
    FILE *fp;

    if (SIM_FLAG == 1) fp = stdout;      // called from command line
    else fp = fopen("results.txt","w");  // called from batch

    fprintf (fp,"\n-----------Simulation results-------------------\n");
    fprintf (fp,"Cutoff=%.0e%s\n",Pz.cutFlt,CpMtype[Cp.useFlag]);
    fprintf (fp,"Contigs= %d\n",n_Zset); 
    fprintf (fp,"F+= %d\n",false_pos); 
    fprintf (fp,"F-= %d\n",false_neg); 
    fprintf (fp,"Gaps= %d\n", num_gaps);

    for (i=0;i<n_Zset;i++) {
        if (Zset[i].n_Qs>0) {
           num_Q += Zset[i].n_Qs;
           num_Q_ctg++;
        }
    }
    fprintf (fp,"Qs= %d, Q contigs= %d\n", num_Q, num_Q_ctg);
    fprintf (fp,"Out of Order= %d\n", compute_out_of_order());
    fprintf (fp,"Chimeric contigs= %d\n",compute_chimerics());
    fprintf (fp,"Singletons= %d\n",count_singletons());
    for (i=0;i<n_Zset;i++)
         avg_score += Zset[i].score;

    fprintf (fp,"Average CBMAP score= %f\n",avg_score/(float)n_Zset);
    fprintf(fp,"---------End of Simulation Results ----------------\n");
}

//This function initialises the required parameters for simulation results
static void init_simulation()
{
   false_pos = 0;
   false_neg = 0;

   if (SHARED_OPTION) {
      sem_init(&false_pos_mutex,0,1);
      sem_init(&false_neg_mutex,0,1);
   } 
   num_gaps = 0;
}

/* This function computes the number of out of order 
   The clones are ordered based on coordinates, then adjacent clones must have
   adjacent names. The one caveat is the if two clones start in the same place,
   their order is indeterminate. This can be fixed in sortctg of clonegen.c,
   before returning zero, check to see if SIM_FLAG and if so, check names.
   Of course, we still would need to know whether the clones are in ascending
   or descending order.
*/

static int compute_out_of_order()
{
   int i=0;
   int num_out_of_order = 0;
   int fp;
   struct contig * cptr;
   CLONE *clp1,*clp2;
   extern int sim_is_adjacent();
   struct remark * tmp;

   for(i=1;i<=n_Zset;i++)
   {
       find_Clone(i);
       orderclones(i);
       for(cptr=root;cptr != NULL && cptr->new!=NULL;cptr=cptr->new)
       {   /*"new" is next; "next" is clone index*/
          clp1=arrp(acedata, cptr->next, CLONE);
          cptr = cptr->new;
          clp2=arrp(acedata, cptr->next, CLONE);

          for (fp=0, tmp=clp1->remark; tmp != NULL && !fp; tmp=tmp->next)
            if (strstr(tmp->message,"FPlus")!=NULL) fp=1;
          if (fp) continue;
          for (fp=0, tmp=clp2->remark; tmp != NULL && !fp; tmp=tmp->next)
            if (strstr(tmp->message,"FPlus")!=NULL) fp=1;
          if (fp) continue;

          if (!sim_is_adjacent(clp1->clone,clp2->clone) && abs(clp1->x - clp2->x)>5)
          {
             //printf("oop %s %s\n", clp1->clone,clp2->clone);
             num_out_of_order ++; 
          }
       }
   }
   return (num_out_of_order);
}
    
/* This function computes the number of chimeric contigs */
static int compute_chimerics()
{
   int i;
   struct contig * cptr;
   CLONE *clp1;
   int found_chimeric = 0;
   int chi_count = 0, stack_count=0, fp; 
   int in_stack=0, left_end;
   struct remark * tmp;

   for(i=1;i<=n_Zset;i++)
   {
       left_end = in_stack =  found_chimeric = 0;
       find_Clone(i);
       orderclones(i);
       for(cptr=root;cptr != NULL ;cptr=cptr->new)
       {   /*"new" is next; "next" is clone i ndex*/
          clp1=arrp(acedata, cptr->next, CLONE);

          for (fp=0, tmp=clp1->fp_remark; tmp != NULL && !fp; tmp=tmp->next)
            if (strncmp(tmp->message,"FPlus",5)==0) fp=found_chimeric=1;

          if (fp) 
          { 
              in_stack=1;
              left_end = clp1->y;
          }
          else {
             if (in_stack) {
                if (abs(left_end - clp1->x)> clp1->fp->b2*2) {
                   in_stack=0;
                   stack_count++;
                }
             }
          }
       }
       if (in_stack) stack_count++;
       if (found_chimeric) {
          chi_count ++;
          // printf("   Ctg%d FP stacks ~%d\n",i, stack_count);
       }
    }   
    return (chi_count);
}
